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A multi-layer model for self-propelled disks 
interacting through alignment and volume exclusion 

Pierre Degond* Laurent Navoret^^ 


Abstract 

We present an individual-based model describing disk-like self-propelled particles 
moving inside parallel planes. The disk directions of motion follow alignment rules 
inside each layer. Additionally, the disks are subject to interactions with those of 
the neighboring layers arising from volume exclusion constraints. These interactions 
affect the disk inclinations with respect to the plane of motion. We formally de¬ 
rive a macroscopic model composed of planar Self-Organized Hydrodynamic (SOH) 
models describing the transport of mass and evolution of mean direction of motion 
of the disks in each plane, supplemented with transport equations for the mean 
disk inclination. These planar models are coupled due to the interactions with 
the neighboring planes. Numerical comparisons between the individual-based and 
macroscopic models are carried out. These models could be applicable, for instance, 
to describe sperm-cell collective dynamics. 


1 Introduction 


Collective motion in systems of self-propelled particles is the subject of a vast literature. 
How collective motion emerges from the underlying local interactions between the agents 
is still poorly understood. The interactions are either of cognitive nature (such as in 
birds, mammals [^) and/or are mediated by a surrounding fluid (such as in swimming 
bacteria, sperm cells, etc. (^). There are several competing strategies to model collective 
dynamics. One strategy relies on individual-based models (e.g. [H[5|-[8, 17,19,21,^) that 


describe how the position and velocity of each individual evolves in the course of time. 


Another strategy relies on continuum models (e.g. [^[3,25,^) describing the system by 
locally averaged quantities such as the mean density, mean velocity, etc. An intermediate 
category of models consist of kinetic models 15 which describe the individual motions 
in a probabilistic way. These different types of models can be connected one to each 
other as kinetic models can be seen as resulting from an inhnite particle number limit 
of individual-based models, while macroscopic models follow from a hydrodynamic or 
diffusive rescaling of the kinetic models and subsequently passing to the limit of a large 
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rescaling factor. We refer to Ref. [1^ for an illnstration of this methodology in the case 
of the Vicsek model (see also below). 

In this paper, we are concerned with hnding a snitable modeling framework for the 
three dimensional motion of spermatozoa in the seminal plasma. As there are evidence 
that sperm-cell motion in the most common experiments is mostly planar 24 , we propose 


a mnlti-layer model where the motion of sperm-cells is planar and sperm-cells may interact 
with other sperm-cells of the same layer or of the two neighboring layers. Sperm-cell 
concentration in raw sperm is incredibly high (as large as 5 10® cm“^) and the nse of 
individnal-based models to reprodnce real sperm-cell experiments is intractable. It is 
therefore necessary to derive a macroscopic model describing the collective motion of the 
cells within each layer and their interactions with the neighboring layers. 

Spermatozoa can be assimilated to two-dimensional discs. Indeed, as regards the occu¬ 
pied volume, flagella can be neglected and spermatozoa reduced to their head. Although 
the heads resemble flat ellipsoids, we simply model them as inhnitely thin flat discs. The 
acting flagellum produces almost constant propulsion. Thus, it is a good approximation 
to suppose that all sperm-cells move with the same constant speed and that only the 
velocity direction is subject to changes. Finally, each disk possesses some inclination with 
respect to its plane of motion, measured by an inclination angle. Therefore, the position, 
velocity and attitude of each disk can be described by the position of its center of mass, 
its velocity direction and its inclination angle. 

Interactions between sperm-cells are mostly hydrodynamic interactions (due to the 
perturbation of the fluid velocity induced by the motion of the cells) and volume exclusion 
(or steric) interactions (due to the impossibility that two sperm-cells overlap). Modeling 
hydrodynamic and steric interactions within dense suspensions of active particles is a 
difficult subject. However, it has been shown in that for self-propelled elongated 
particles, such interactions simply result in local alignment of the particles with their 
neighbors. Therefore, we assume that all these interactions can be lumped into a local 
alignment interaction with the close neighbors. 

As already mentioned, we assume that particle motion is planar and takes place in 
parallel two-dimensional layers. Each particle belongs to one layer for all times without 
the possibility to change its layer. For the reasons outlined above, spermatozoa interact 
inside these 2D layers by alignment of both their velocity and inclination with those of 
their close neighbors. Specifically, we consider the time-continuous version of the Vicsek 
microscopic model as proposed in Ref. 12 , where each particle tends to align with the 
mean direction of its neighbors up to a small Brownian perturbation. The original model 
was proposed by Vicsek et ah 28 and several variants have been proposed in |^pT|[l6|[2^ . 


The inclination alignment dynamics follows a similar rule with the exception that the 
interaction is nematic (i.e. two inclination angles differing by a multiple of vr lead to 
the same disk attitude). The combined alignment dynamics in the velocity-inclination 
variables thus differs from the 3D Vicsek dynamics (specifically, in the 3D Vicsek, velocity 
belongs to a 2D-sphere whereas here the pair (velocity, inclination) belongs to a 2D torus). 

The different layers also interact via the volume-exclusion constraint. Indeed, due the 
inclination of the disks, the spermatozoa of one layer exert a friction on the spermatozoa 
of the nearby layers. This interaction results in increasing or decreasing the inclination 
of the spermatozoa in these layers. A given layer thus acts on the neighboring ones in 
a similar way as the wind does on plant canopies (Til. The involved mechanical forces 
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depend on the geometric confignration of the discs: it thns depends on their respective 
inclinations and also on their respective velocities. It results in alignment of velocities 
and repulsion of inclination angles. Layers are thus coupled and this coupling depends on 
the so-called overlap function that quantihes the distance between layers. 

We then consider a mean-held kinetic version of the model. This equation provides 
the time evolution of the distribution function in phase space (position, velocity and 
inclination) and takes the form of a Fokker-Planck equation. Here, it is formally derived 
in the limit of an inhnite number of interacting particles. Although the mathematical 
validity of the mean-held limit has been established for the Vicsek model |^, it is still 
open for the present model and our result so far is only formal. We perform a spatio- 
temporal hydrodynamic rescaling of the mean-held kinetic model, considering that the 
intra-layer interaction scales are much smaller than those of the inter-layer interactions 
and that the latter occur at the same scales as the macroscopic evolution of the system. 
There results a singularly perturbed Fokker-Planck equation, involving a small parameter 
£ measuring the ratio of the small (microscopic) scale to the large (macroscopic) one. 

The macroscopic description of the system is found by letting e: to zero in the singularly 
perturbed mean-held kinetic model. We hrst need to hnd the equilibria associated to 
the Fokker-Planck operator. We show that these are given by products of von Mises 
distributions in the velocity and inclination angles respectively, von Mises distributions 
are the natural analog of Gaussian distributions for probabilities on the sphere. We then 
need to integrate the equation against the collisional invariants. However, as noticed in 
Ref. 12 , only mass is a collisional invariant and we are thus lacking two more collisional 
invariants to obtain the dynamics on the velocity and inclination. Following Ref. 12 


we have to introduce the “generalized collisional invariants” of this operator: these are 
collisional invariants valid only on functions with prescribed mean velocity and mean 
inclination. We then are able to derive the macroscopic model. 

The obtained model consists of a continuity equation for the density p and two evolu¬ 
tion equations for the mean velocity angle p G [0, 27r] and mean inclination angle 9 G [0, tt]. 
All these quantities are indexed by h corresponding to the layer. The model is written: 


dtPh + CCi V:,; • {phV{ph)) = 0, 

Ph{dt^h + cc2{V{iph) ■ '^x)Ph) H- y{^h)^ ■ ^xPh 

Kl 


C3 


{gM2M2){9,9k)c^pkV{pk). 

fc, k—h=^l 

phidtOh + cci{V{(fih) ■ '^x)9h) 


( 1 ) 


( 2 ) 


C 4 


{cipV{iph))^ ■ ^ sgn{k - h){gM2M2d0l2){O,9k)cipkV{ipk), (3) 


k, k—h=±l 


where V{(p) = (cos sin (^)^ denotes the velocity vector. The constants will be dehned 
further. The left-hand sides of ([^-([^ form the SOH (Self-Organized Hydrodynamics) 
model describing the Vicsek dynamics at the macroscopic level 


12 


They respectively 

account for the conservation of mass and convection of the mean velocity angle. The 
left-hand side of ([^ describes the advection of the inclination with the same advection 
velocity as for the mass. Finally, the right-hand sides of ([^-([^ describe the inter-layer 
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interactions. The obtained model is thus an hyperbolic system (like the SOH model) with 
source terms that couple the velocity and inclination dynamics. 

We remark that, once the velocities of all the different layers are co-linear, the source 
terms vanish. The model thus simplihes into a superposition of standard SOH models 
in each layer. However, before reaching an equilibrium, the interplay between inclination 
and velocity crucially determines which equilibrium velocities and inclinations will be 
attained. 

To validate the macroscopic model, numerical simulations are performed and compared 
with those of the particle model. With this aim, we adapt the numerical relaxation 
method of 20 designed for the Vicsek model. The numerical simulations show that the 


macroscopic model captures the velocity alignment between the different layers quite well. 
However, some differences in the inclination dynamics appear. Indeed, some transient 
“meta-stable” conhgurations arise during the course of time. They are more rapidly left 
away by the microscopic dynamics than by the macroscopic ones, probably because of 
the stochastic fluctuations associated with the microscopic dynamics. This effect due 
to the hniteness of the particle number in the microscopic dynamics could probably be 
reproduced by including a stochastic term in the macroscopic model such as the one 
proposed in Ref. [2^ . 

The outline of this article is as follows. In Section we introduce the microscopic 
model. In Section we present the mean-held limit and the hydrodynamic scaling. In 
Section we provide the derivation of the macroscopic model. In Section we compare 
the microscopic and macroscopic dynamics on several test-cases. A discussion of the 
results is provided in Section]^ [A| provides complements on nematic alignment modeling, 
[B] gives some properties of the coefficients of the macroscopic model and describes the 
numerical schemes used for the simulations. 


2 Microscopic model 

Spermatozoa are represented by discs of radius R moving in different layers. In hrst 
approximation, the layers, indexed by Z, are copies of the plane: layer h denotes the 
plane {(x, ?/, z) G z = hd] where d > 0 denotes the inter-distance between the layers. 
With no interactions between layers, spermatozoa are supposed to be orthogonal to the 
layer and follow the Vicsek dynamics. However, if the layer inter-distance is lower than 
the disc radius R, spermatozoa of layer h will exert a force on the spermatozoa of the 
neighboring layers h — 1 and h + 1: they may force them to incline (with respect to the 
layer plane). 

We consider N discs in labeled by k G {1,...,V}: each disc is contained into one 
layer and thus disc k has a permanent altitude hk G Z. The two-dimensional movement 
into the layer is described by the position of its center of mass Xk{t) G and the velocity 
orientation of its center of mass 14(f) G we indeed consider that all particles move at 
the same speed c > 0. The velocity of the particle is thus given by cI4. We introduce 
the angle (pk(t) of 14(f) with respect to a reference axis, so that 14(f) = V{ipk(t)) with 
V((p) = (cos (p, sin (p). 

Concerning the conhguration of the disc in space, we suppose that the disc moves in 
one direction contained in its plane: the orientation 14 (f) ^ belongs to the disc plane. 
The angle of this plane with respect to the z-axis is denoted 0k{t) G ]R/[0,7r]. An angle 
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(b) One spermatozoon 


Figure 1: Spermatozoa in layers 

6k = 0 or 71 means that the disc is perpendicular to the plane while an angle 6k = ±7i/2 
means that the whole disc lies in the layer plane. 


2.1 Dynamics for the centers of masses. 


The centers of masses follow a Vicsek-like dynamics as introduced in Ref. 12 . The 


dynamics of the positions and the velocities are given by the following equations: 
dXk 


dt 


-{t) = cV{ipk{t)), 


d(pk{t) = sm{ipk{t) — (p^k^{t))dt + a/ 2D dBf’^ modulo 27r. 


( 4 ) 

( 5 ) 


Two dynamics are in competition: alignment and diffusion. Each particle of a given layer 
hk tends to align with a direction with an intensity z/ supposed constant. The 

direction V{ip)^^{t)) is dehned as a weighted mean direction of the neighbors particle k in 
layers hk, ± 1 within the disc of radius Ri: 






m = E 

j,hj=hk,\Xj{t)-X,,(t)\!^Rx 


^V^,w,nb(^) ^ E 9{.0j{t),6k{t))V{ipj{t)), 

y =h.j. it 1,1 Xj (P — JVfe (t) I sgRi 


( 6 ) 

(7) 


where (f) denotes the contribution of neighbors belonging to the same layer hk, 
denotes the contribution of neighbors belonging to layers ± 1 and f3 quantihes their rel¬ 
ative involvement. Superscripts “nb” means neighboring layers and “w” means weighted. 
Indeed, due to steric constraints within layers, directions of neighbors of layers ± 1 are 
weighted according to their inclination. Supposing distance h between layers is smaller 
than twice the particle radius, h < 2R, we dehne the overlap function g by: 

ff(dj,6k) = ^ (R(| cosgji + |cos4|) -^) + - (8) 

This function quantihes the overlapping area of two discs in the 5 direction. The particle 
direction is also submitted to a Brownian motion with diffusion coefficient D > 0. 
We here neglect congestion forces and we also neglect alignment between discs of different 
layers. 
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2.2 Dynamics of the disk orientations. 


The disc angle dynamics follows the following torqne balance equation in the over-damped 
regim^ 


dOkit) = {—K sin(2(6lfc(f) — Ok(t))) + Tk(t))dt -|- dB^’^ modulo tt, (9) 


where the alignment dynamics inside each layer is in competition with steric forces be¬ 
tween layers. The factor 2 inside the sine function takes into account that the inclination 
interaction between disks is a nematic one, i.e. orientations 6k = 9k or 9k = —9k are 
equivalent. In this way, Eq. ([^ preserves the fact that 9k is dehned modulo vr. K is the 
rotational stiffness and 5 > 0 is a diffusion coefficient. 

The inclination of each particle in layer hk tends to align with the mean inclination 
9k of neighboring particles of the same layer in the disc of radius R 2 . The mean inclina¬ 
tion is dehned through the following average, which corresponds to a nematic alignment 
mechanism (see 0: 

4 m- Y. (“) 

' j,hj=h^„\X,{t)-Xk{t)\^R2 

where 9k(t) is dehned modulo vr. 

Particles on neighboring layers ± 1 exert a steric force on particles of layer hk- Tk{t) 
is sum of the weighted torques Tkjit) (with respect to the longitudinal axis) of the forces 
exerted by the discs: 

Tk{t)= Y. g{9,{t),9k{t))Tkj{t), ( 11 ) 

j, hj=hk±l, |Xj(t)-Xfe(t)|sgi?3 

Tkj{t) = /ii?sgn(hj - hj) iVkit) x Vj(t)) ■ z, (12) 

where ;§ is the unit vector in the vertical direction and g{6j,9k) is the weight dehned in 
(j^. Indeed, the force Fkj exerted by disc j on disc k is supposed to be the projection of 
the velocity direction Vj on the orthogonal plane to 14 : 

Fkjit) = 

where /r is a mobility coefficient. Then the torque Tkj of the force Fkj with respect to the 
longitudinal axis of the disc is given by: 


n,m = i(Xj(i) - x,(i)) X n,(i)l • mm 

and using the approximation {Xk{t) — Xj{f)) ~ Rsgn{hk — hj)z, we obtain: 

Tkjit) ^ laRsgnihk-hj)[z X Py±(^t)Vjit)]-Vkit) 

^ jiPsgnihk - hj) [Py±(i)V,(f) X Vkit)] ■ z 
^ /iPsgn(hfc - hj) [Vjit) X Vkit)] ■ z 

^The torque balance equation reads: I fjft + -I- K sin(2(0fe — 4)) =Tk, where I is the moment of 

inertia of the disc with respect to its longitudinal axis (parallel to 14) and C is the dissipation coefficient. 
In the over-damped regime, I is negligible and we recover ([^ with K = K/C. 
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3 Mean-field kinetic model and rescaling 

3.1 Mean-field kinetic model 

We introduce the distribution function in phase space: f{x,(p,9,h,t), where / is 27r- 
periodic in ip and tr-periodic in 6. Note that h is still a discrete parameter numbering the 
layers. The distribution function satishes the mean-held kinetic model; 

a,f + cv, ■ (V(<p)f) = sm(2(« - },)) + T,)f) + sa}f 

-a,(-i^sin(^-45}“)/)+£'S5/. 

where 9f = 9f{x, h, t) and = (py^{x, h, t) are dehned by: 




jR^jix,h,t) = 


/ 0 G[O, 7 r],(^G[O, 27 r], 

yeM2,|j;-x|^R2 


e * fiy,p,9,h,t)dydpd9, 


and 


V{pf{x,9M) = 


Kf{x,h,t)+l3J‘fl'”f'°{x,9,h,t) 


\Jr^j{x, h, t) + f3J‘^’^f°{x, 9,h,t)\' 


J‘^^j{x,h,t) = 


/ 0 G[O, 7 r],i^G[O, 27 r], 
y&R'^ ,\y-x\^Ri 


V V?, h, t) dy dp d9, 


k—h=:kl 




fO^ [0,7r] [0,27r], 

y&R‘^,\y-x\^Ri 


g{9, 9')V{p)f{y, p, 9\ h, t) dy dp d9'. 


The torque Tj = Tf{x,h,t) is given by: 

Tf{x, p, 9, h, t) = [Nf{x, 9, h,t)xV (y?)] • 

Nf{x,9,h,t) = pR sgn{k — h) J‘^^jr{x,9,k,t), 

k, k—h=:kl 

and depends on neighboring layers h — 1 and h + 1. 


3.2 Hydrodynamic rescaling 

We then perform a hydrodynamic rescaling to look at the large time and space scale 
dynamics. The hydrodynamic rescaling consists in introducing macroscopic variables in 
space and time: x' = ex, t' = et, with e 1. After dropping the tildes, the kinetic 
distribution f^{x', v, w, h, t') = f{x, v, w, h, t) satishes the following equation: 

e{dtr + cV. . (r(v>)/')) = -S9((-A'siri(2(» - 0}.)) + T^f) + Sd^f 

+ ap,^smCp-vT’)r)+Dalr, 


(13) 
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where the rescaled Oj{x, h, t), h, t), Tj{x, h, t) are dehned by: 




^ ’ ’ ’ \JeRU^^^h,t) + /3J^^^f{x,e,h,t)\' 

Tf{x,e,ip,h,t) = [Nj{x,e,h,t) X V{lp)\ ■ z, 

AfJ(a;,6',h,t) = ^ sgn{k - h) J^£j{x,9,k,t). 

/c, k—h=:kl 

We can easily show the following expansion: 

JsRj = e‘^{nR^jf + 0(£^)), 

with 

jj {x, h,t) = f V((p)f(x, ip, e, h, t) dp d9, 

J 0G[O,7r], £/?G[0,27r] 

9, k,t)= [ g{9,9')V{p)f{x, p, 9', k, t) dp d9', 

J 6S[0,7r], (/jS[0,27r] 

jj(x, h,t) = f e^'-^fix, p, 9, h, t) dp d9, 

J 0G[O,7r], ;/:G[0,27r] 


(14) 

(15) 

(16) 
(17) 


(18) 

(19) 

( 20 ) 


where the last qnantities are the localized mean inclination angle and momentnm. There¬ 


fore, system (13 17) becomes, dropping 0{e‘^) terms: 


where 


e{dtr + cV. • iVip)n + deiTf.n) =Kdg{smi 2{9 - 9 f.))r) + 

+ izd^{s\n{p-p^^^)r)+Ddlr, 


2r5pxXt) ^ 


( 21 ) 


( 22 ) 


\i'f (x,h,t)+13 ij'”'"'’)!, 0,h,t)\' 
fp-'''’{x,0,h,t) ^ 5^ jf'’(x,e,k,t), 

k, k—h=:kl 

Tf{x, p, 9, h, t) = [Nf{x, 9, h,t)xV ((p)] • z, 

Nf{x,9,h,t) = - pRnRl sgn{k — h) jj''”{x,9,k,t). 

k, k—h=^l 


(23) 

(24) 


We snppose also that the interaction between the layers happens on large time scale. 
Therefore, we write: ^pRttR\ = p' = 0(1). Inserting this ansatz, we obtain: 

Nf{x,9,h,t) = p' ^ sgn{k - h)j‘^’'"{x,9,k,t). 

k, k — h=:kl 


(25) 
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Moreover, we suppose that (3 = e(3' with (3' = 0(1). We thus have the following expan- 
sioi£l 

I ^,w,nb/ 0 f)\ 

+ £/?' 0, h, t)) + 0{e\ 


where 93 / and are dehned by: 

V(- ( h T^r-w,nb. p , 

V{(pf{x,h,t)) = ^ V{(pf’ {x,e,h,t)) = 


\jf{x,h,t)\ 


\jf-’-\x,e,h,t)\^ 


(26) 


and Pv[^f{x,h,t))^X denotes the projection of X onto the orthogonal plane to V{(pf{x, h, t)). 
Taking the cross product of the previous expansion with V{ip), we easily obtain: 




(/?,w,nb I 
/' 


sin(93 - ipj^) = sin(93 - ipfe) - e(3' sin((y9p“ - (pfe) cos{(p - (pf.) + 0(£^). 

\J f^\ 


Therefore, regrouping the 0{e) terms together, equation (21) becomes: 

e(d,f + cv, ■ (r(^)/') + a,(T,.r) + aps,.r)) 

= K 59(sm(2(« -»/.))/') +Sdlf + vdPsm(ip-tpi,)S’) + Dd^f, (27) 


with 


\3. 


(/?,w,nb I 
/ 


Sf. = iyl3' '‘'f ' sin((p7f’^ - cos{p> - ipp). 

\3p\ 


(28) 


System (22)-(28) is the starting point for the derivation of the macroscopic model. 


4 Macroscopic model 


4.1 Equilibria 


We want to take the limit e: —)• 0 in system (22)-(28). Therefore, assuming that the 
distribution function converge to a limit denoted by as e ^ 0, this limit satishes 
the equilibrium condition Q{f^) = 0 where 

Qif) = K dg{sm{2{e -Of)) f) +6d‘^f + iyd^{sm{ip-ip f)f) +Ddlf. 


^For any vectors a, 6 S and e > 0, we have : 
a -I- e6 a -I- e6 


a 


+ sb\ 


+ £^-5 + 0 ( s 2 ) 


• |a| 

£ ( Id- 


a 0$ a 


= (o + £o)(i-^o-o + o(£)) 


+ 0{e^). 
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We define the von Mises-Fisher (VMF) distribntion with periodicity 27 i/n, n G N\{0}, 
concentration parameter /t > 0 and direction -00 £ M/(27rnZ) by: 

= ^^exp (Kcos(n('0 - V’o))), (29) 


with 


We introduce 


Zn,K = / exp (fi;cos(n(V’ - V’o))) dip. 

Jo 




(30) 


with 


z/ 




We show the following: 

Proposition 1. We have 


Qif) = 6de[M^.eMe)de 






+ Dd^(ip, 9) 0^))- (31) 


Proof. We have 


Qif) = de{K sm{2{e - 9 f))f + 6 dgf) + d^{iysm{(p - (ff)f + D d^f) 
= ddg{- ^dg{cos{2{9 - 9f)))f + dgf) 

+ Dd^{- ^dp!{cos{p! - ^f))f + d^f). 


In view of (29), Eq. (32) can be written: 

f{p,o) 
^2^,0,id) 


0 = 6dg(M,K.{9)dg 


+ Dd^[Mi^^^^^{(p)d^ 


f(p,o) 


and from (30), we deduce (31). 
We dehne 

Pf{x, h,t) = 


' 6S[0,7r], (/jE[0,27r] 


/(x, ip, 9, h, t) dp d9. 


(32) 


(33) 

□ 


We have the: 
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Proposition 2. For any distribution we have 

u K - 

P[pM^,o] = = Ci(l, n) 3[pM,^,] = Cl(2, -j) 


with 


In particular, we have: 


ci(n,/t) = 


D 


codnu) du. 


JO 


^[pM^g] = P modulo TT, PipM^g] = modulo 27 r. 


(34) 

(35) 


This proposition shows that for the distribution pAi^^g, p is its local density and 99 , 6 
are its mean direction of motion and mean inclination respectively. 

Proof. By easy computations, we have 

^[pM-g] = [ 0) dip de 

0G[O,7r], 99G[0,27r] 

= p(/ M,„(B)de)(f V{<p)M,,~,f(,p)d<p) 

^J9£[0,n] ^ ^ ^ J</je[0,27r] ' 

= pV{(p) / Mi^j^^^{ip) COS (p dip 

J (pS[0,27r] 

= Cl(l,^)pP(<^), 

We note that ci(n, k) is such that 0 < ci(n, /t) < 1. Therefore, 

V(-f\pM,P = = V(v), 

'd[pM.^ g\ I 

showing that PipM^g] = P rnodulo 27r. Similarly, we have 

dipM^^g] = ci(2, 

and so showing that 0\pM^ e\~d modulo tt. □ 


Proposition 3. The set £ of equilibria, i.e. £ = {/((p,^) | / >0, and Q{f) = 0} is given 
by : 

£ = [pAi^g I (p e [0,27r], P e [0,7r], p e M+} . 
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Proof. From (31), we deduce that 

/ 


/ M 

'0e[O,7r],(,ae[O,27r] 


iff,Of 


dip dO 


d 


de{ 


f 




PD 


J 0e[O,7r], (/je[0,27r] 

Now, let / be an equilibrium, i.e. Q{f) = 0. Then, 


d^{ 


f 




^ 6S[0,7r], iy3S[0,27r] 


+D\d^{ 




f 




) ) M^^^g^dpd9. 

(36) 

ppPfdp d9 = 0. 

(37) 


This shows that there exists p G M+, such that / = which implies that there 

exist p G M+, G M (modulo 27r) and 0 G M (modulo vr) such that / = pM.cp,e- 
Reciprocally, we show that / = pM.cp,e s-re such that Q{f) = 0. Indeed, this follows from 
( [^ and (37) if we show that for / = pAi^pp, dj = 9 and pf = p. But this follows in turn 
from Proposition □ 

According to this proposition, there exists p{x,h,t) G and p{x,h,t) G [0,27r], 
9{x,h,t) G [0,7r] such that: 

lim/^ = /°, with f{x,p,9,h,t)=p{x,h,t)Mp^^^h,t)Ax,hpiP^^)- (38) 

Now, the goal is to hnd equations for (p, <p, 9) as functions of {x, h, t). 

4.2 Collisional invariants 


Eq. (27) can be written 


Si/' + cV. . (l/(v>)/') + S»(T/./') + 3„(S/./')) = -Q(/'). 


(39) 


In the limit £ —)■ 0, the right-hand side of (39) is singular. The goal of a Collision 


Invariant is to cancel this singular term through integration in (<p, 9) against a suitable 
test functions. For this purpose, we dehne: 

Definition 1. A Collision Invariant (Cl) is a function I{p, 9) such that for all function f{p, 9), 
we have 

[ Q{f)Idpd9 = 0. 

J i/?G[0,27r] 

We denote by C the space of Cl. It is a vector space. 


Here, clearly, C contains the constant functions, since, as seen from (31), we have 


/ Q{f)dpd9 = 0. However, no other Cl appears obviously from (31). The use of the 


constant functions as Cl already gives the mass conservation equation. Indeed, integrating 


(39) with respect to {p,9), we get 


dtPfs + Vx • jje — 0. 


(40) 


In this equation, the 1/e singularity has disappeared and the limit £ —)■ 0 of Eq. (40) leads 


to an equation for /°. However, as seen from (38), /° depends on three scalar quantities 


p, 9 and p but (40) is only one single scalar equation. Therefore, it is not sufficient to 
determine the dynamics of /°. For this reason, we look for a weaker invariant concept, 
that of Generalized Collision Invariant, as defined in the next section. 
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4.3 Generalized Collisional Invariants. 

For a given pair {(p,9) G M./2nZ x M/ttZ, we define the operator Q{(f,9;f) by 

^na,{M,,ir,B)a,(J^)). (4i) 

We note that / h-)■ Q((f, 9; f) is a linear operator and that 

Q{f) = Q{^f,9f,f). (42) 

We define the Generalized Collisional Invariants by the following 

Definition 2. Let {cp, 9) E ]R/ 27 rZ x M/ttZ be given. A Generalized Collision Invariant (GCI) 
associated with {p,9) is a function I{p,9) such that 


Ge[0,7r], Iy3e[0,27r] 


Q{(p,9-, f)Idipd9 = 0, 


TT 


V/ such that 9f = 9 mod. — and ipf = ip mod. vr. 
We denote by Q the space of GCI associated with {p,9). It is a vector space. 


Referring to (29), for the simplicity of notation, we define 

^2,e = ^2,K2,e- 

We introduce the two following functions: 

(i) The function Ii{p) is a 27r-periodic solution of the problem 


(43) 


(44) 


= sinpMi^o, / Ii{p) dp = 0 

J t/3S[0,27r] 

(ii) The function l2{9) is a vr-periodic solution of the problem 


(45) 


ae(kh„ ash) = siii(29)M2,„, / hW d0 = 0. 

J 0G[O,7r] 


(46) 


Now we have 

Theorem 4. The solutions R and I 2 are unique. Moreover, the space Q of GCI associated 
to {Pt9) is three dimensional and spanned by 

X,{p,9) = l, X^{p,9) = h{p-p), X2{p,9) = l2{9-9). 
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Proof. Introducing the L^-adjoint 9] /) of Q{(p, 9] /), we can write: 

[ Q{if,9-f)Xdipd9= [ fQ*{ip,9-X)difd9. 

J 0G[O,7r], (/:G[0,27r] J 0G[O,7r], (/?G[0,27r] 

The constraints 9f = 9 modulo | and cpf = cp modulo vr can be equivalently written: 
f f sm{2(9 — 9)) d9 dp = 0, f f sm{p — p) d9 dp = 0. 

J 6s[0,7r], 93g[0,27r] J 0g[O,7r], ipS[0,27r] 


Since these are linear constraints, by a classical duality argument, (43) is equivalent to 
saying that there exist (/3‘^,/9®) G such that 


^0S[O,7r], yJS[0,27r] 


/ Q*{p,9]X) dpd9 


f (/3^ sin(2(6' — 9)) + (3'^ sin((y9 — p)') d9dp. 


>^0S[O,7r], 95S[0,27r] 

for all functions / without constraints. This implies that X satishes: 

G such that Q*{p,9]X)=(3^ sm{2{9-9))+/3'^ sm{p-p). 

or, using the explicit expression of Q*-. 

5 de {M^,e del) + D {M^^g d^X) 

= (/3® sin(2(0 - 9)) + [3^ sin((p - p)) M^^g. 


(47) 


Multiplying by a test function J', integrating with respect to {p,9), using Green’s for¬ 
mula and the 27r periodicity in p (resp. the vr periodicity in 9), we obtain the following 
variational formulation: 

[ M^pg{ 6 dgXdeJ + Dd^Xd^j)dpd9= f 'ipjdpd9, (48) 

J 0S[O,7r], 93S[0,27r] J 0S[O,7r], 93S[0,27r] 

with xf = — (/9® sin(2(6' — 9)) + (3"^ sin((p — p))XA^^ 0 . We now introduce the functional 
spaces L^([0,27r] x [0,7r]), i7^([0,27r] x [0,7r]) endowed with their classical Hilbert norms 
and inner products, together with: 

Hi ([0,27r] X [0,7r]) = {J e H\[0,27 i] x [0,7r]) | J{0,9) = J7(27r,0), 


J{p,0) = J{p,7r), a.e. {p,9) G [0, 27r] x [0,7r]}, 


and 

Hi([0,27r] X [0,7r]) = {H;,,([0, 27r] x [0,7r]) | f Jdpd9 = 0}. 

J 0S[O,7r], y3S[0,27r] 

Thanks to a PoincarGWirtinger inequality (which can be easily proved using the Rellich- 
Kondrachov compactness theorem), the semi-norm 
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is a norm on 2tt] x [0, vr]) equivalent to the classical norm. Therefore, thanks to 

Lax-Milgram theorem in hf^([0, 2tt\ x [0, tt]), there exists a unique X G hf^([0, 2tt\ x [0, tt]) 
such that (48) holds for any J G if^([0,27r] x [0,7r]). Now, since ip is the sum of two 
terms, one being odd in 6 — 6, the other one being odd in ip — cp, we have 


' 0E[O,7r], (^G[0,27r] 


Ip dpdO = 0. 


Therefore, X satishes (48) for all J G H^^^{[0,27r] x [0,7r]) (and not only for J G 


hr^([0,27r] X [0,7r])). Furthermore, all solutions in H^^^{[0,27i] x [0,7r]) of (48) equal the 


unique solution in hf^([0,27r] x [0, tt]) up to a constant. Indeed, if X solves (47) with -0 = 0 
in i7per([0,27r] X [0,^, we have \X\\ = 0 and therefore, X is a constant. 

Now, we solve (47) for = {D,0) or (0,5) and for this purpose, we use the 


functions Ji and I 2 dehned at (45) and (46). We note that = h{x ~ x) is the 

unique solution in i7^([0,27r]) of the variational formulation 

I d}(pl\^(p dp 

J ip£[0,2n] 

= — I sin{p — p) Jidp, VJi G i7^([0, 27 r]). (49) 

J t/3S[0,27r] 

and that l 2 p{d) = l 2 {d — 0) is the unique solution in i7^([0,7r]) of the variational formu¬ 
lation 


' 0 G[O, 7 r] 


^2,e 90^2,9 90 X 2 dO 


' 0 G[O, 7 r] 


sin(2(» - «)) Af 2 ,« J 2 id, VJ 2 6 //‘([O, ir|). 


(60) 


The existence and uniqueness of solutions to (49) and (50) follow from the same kind of 
arguments as for problem (48). Now, it is an easy matter to check that both Ii^cp and 


I 2 P are solutions of (48) with (/?‘^,/3®) = (X),0) and (0,5) respectively. Moreover, they 
both are in 77^([0, 27r] x [0, vr]) and by the uniqueness of the solution of (48), they are the 
unique solution of this problem with these choices of (/d‘^,/3®). We deduce that the space 
Q is three-dimensional, spanned by Xo((p, 6 ) = 1, Xi((p, 6 ) = Ii,,p{p) and X 2 ((p, 6 ) = l 2 p{d), 
which ends the proof. □ 


Functions Ii{p) and hid) equal: 


/.(^) = _- to 

Kl 


du 


hid) = 


2K2 2K2 f 

Jc 

Thanks to (42) and (43), they satisfy: 


d ^tt/2 


cos 2 «^y 


^0S[O,7r], (/JS[0,27r] 


' 0S[O,7r], (/3S[0,27r] 


Q(f) h(v - 'ff)difde = 0, 
Q(nh(d-Sj,)Apde = (i. 


(51) 

(52) 

(53) 

(54) 
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4.4 Macroscopic equations 


We obtain the macroscopic dynamics by integrating system ( 21 ) against the GCI. The 
resulting equations are presented in the following proposition: 

Theorem 5. The density p{x, h, t), the angle of the mean direction Lp{x, h, t) and the mean 
inclination angle 9{x,h,t) satisfy the following system: 


dtp -f cci Va; • {pV(cp)) = 0, 

p(dtp + cc2(V(p) ■ V^)(p} + —V((p)-^ ■ Vxp 


Ki 


C3 


Y. {9M2M2){e,h)ciPkV{pk). 


/c, k—h=:kl 

p{dte + cc^{y{ip)-w^)e) 

p' 


C4 


{cipv{(p))^- Y Sgn{k - h){gM2M2d0l2){O,ek)cipkV{(pk), 


(55) 


(56) 


(57) 


k, k—h=±l 


with 


{gM2M2)il h) = 9(0 + e, h + B')M2{B)M2{B') dBdB', 


(58) 


{3M,_Ahaeh){9,9k) = 2K2 / g(e + OAa 0')AW)deh{ff)M2(e')dede' , (59) 

ie,0'e[O,7r] 

where {pk,Pk,dk){x,t) denotes {p,(p,B){x,k,t) and ci = Ci(1,Ki) (with ci{n, n) defined at 


(35)) and C 2 , C 3 , C 4 given by: 

4g[o,2H sin if COS if Mih dp 


C2 = 


C 3 = 




Gs[0,27r] 


sin pMiIi dip, C 4 = —Ak\ 


^ 6S[0,7r] 


sin 26*M 2/2 dB, 


and Ml = Mi^o. ^2 = ^ 2,0 (see ( 0 ). 

Before giving the proof, let us make some comments. The left-hand side of equations 


(55)-(56) is exactly the SOH (Self-Organized Hydrodynamics) model describing the Vic- 
sek dynamics at the macroscopic level (see Ref. [^). The right-hand side of equation 
describes the alignment of V{p) toward a linear combination of the velocities of the neigh¬ 
boring layers. The weights of this linear combination depends on the inclination of the 


different layers. Equation (56) describes the advection of the inclination with the same 


advection velocity as for the mass. The right-hand side of eq. (56) finally also evaluates 
weighted alignment terms between layers. 


The weights are given by (58 )-(59). They are integral operators, quadratic with respect 


to the macroscopic equilibria in inclination variable. They are scaled in such a way to be 
bounded quantities (see B- Weights ( [ 5 ^ involve the second generalized invariant. 

Finally, this macroscopic model depends on several coefficients, ci, C 2 , C 3 , C 4 , that are 
all positive and bounded by 1 (see 0. They are all averages of the von Mises equilibria 
(either in velocity or inclination variables) against the collisional invariants. 
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Proof. We apply the moment method: hrst integrate the equation against the collisional 
invariants and then taking the limit £ —)■ 0. 

Mass conservation equation. Here it is just a matter of passing to the limit e —)■ 0 


in (40), using (38), (34), (35). We immediately get (55) 


Velocity and inclination equation. We multiply (27) by If := Ii{(p — cpfs) or 
If := 12(0 — 9fe), integrate with respect to {(p,9) and use (53) or (54). We get 


dtf + cViif) ■ V.r + deiTf.f) + ) Ild^d9 

'0e[O,7r],(,5e[O,27r] 

= 0, fc = l,2. 

In the limit £ —)■ 0, we have (<^/£, 9fe) —>■ {(f, 9) and consequently If —)■ Ik, where Ik stands 
for the same quantities with replaced by pAi^ Q. Therefore, we get 


' 0G[O,7r], £/?G[0,27r] 


9t{pAi<p,e) + cV{(p) ■ x{,pM.cp,e) + d0{TpM-gpM.^^g) + 


d^{SpM^ gP-^^fi)) h d^pd9 = 0 , k = 1,2, 
Tedious but easy algebra leads to: 


( 60 ) 


+ cV(ip) ■ VJpM,f_s) = Mff ^ T,j, 

*J = 1 

with Tij even with respect to if f is even and even with respect io 9 — 9 ii j is even, 

odd otherwise. This gives: 

2cK 

Til = P^- sin(2(6' - 9)) sin(93 - (p)V{(p)-^ ■ VA 
0 

Ti 2 = p—si\i{(p - (p)dtp + csin(v3 - (p)V{(p)-^ ■ V^P 

cu 

+p— sin((y9 - (p) cos((p - (p)V((p) ■ VxP 
2K - - 2cK 

T21 = p—^ sin(2(0 - 9))dt9 -L p—^ sin(2(0 - 9)) cos((p - (p)H((p) ■ Tx9 
0 0 

Qjy 

T22 = dtp + ccos((p - p>)V{(p) ■ VxP + p—sin^((p - (p)V((p)-^ • VxP 


Since Aippp is even in 6^ — 6^ and p — p while Ji is odd in p — p and I 2 is odd m. 9 — 9, we 
get: 


'0e[O,7r], (^g[0,27r] 


dt{pM^^s) + cV(p) ■ Vx(pM^^ 0 ) h dip d9 


6e[0,7r], (/je[0,27r] 

u 


M^gTuIidip d9 


= P 


D 


i sin pMi Ji dip ) dtp 

¥’e[0,27r] 

CP 

+ W ^ , 

+c( / sinpMiJidp) V((p)'‘'■ Vxp. 

2(^e[0,27r] 


sin p cos pMi/i dp ) V (p) • V^p 


(61) 
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while 


^6S[0,7r], 95S[0,27r] 


dtipM^^e) + cVicp) ■ V^ipM^^g) h dip dO 


'6e[0,7r], (^e[0,27r] 

f2A: / 


-^(3,0 ^21 12 dip dO 


= pi —— ( / sin 26 M 2 I 2 doi dtO 

^ 0 '^./flcrnTri ^ 


-T( 


'0s [0,71 

/ sin26'M2/2 r / cos pMi dpi V{p) ■ Vxd 

0S[O,7r] ' 0 ipS[0,27r] ^ 


(62) 


We now treat the last two terms of (60). Using integration by parts we have 


' 0S[O,7r], 93S[0,27r] 


de{TpM^-epM^^e) h dp dO 


' 0S[O,7r], 93S[0,27r] 


0S[O,7r], 93S[0,27r] 

d^{SpM^ epM^fi) h dp dO 


'^pM^ sP-^‘P,e dp dO, 


^ 0S[O,7r], y3S[0,27r] 


SpM^^gP-^pP dp,Ik dp dO. 


Since Ji does not depend on 9 , the contribution of the third term of (60) for k = 1 


vanishes. Let us write (28) as follows: 

SpM^-, = I'ld' | I ifpMl) COs{p - PpM^^g)- 

\dpM^-g\ ^ ^ 


From (19) and and using (34), we get: 




' 0 G[O ,71 


g{0,S + g')M2d0') CtpV{f>), 


(63) 


and then: 

^pM^^s = ^Id' cos{p -p)V{p)^ ■ ^ 


g{e,ek + e')M2de‘ 


A cipkV{pk) 


Cip 


k,k-h=±l ^“'®e[ 0 , 7 r] 

Therefore, after integration by parts, and using that Mid^Ii is even in (p, we have: 


^SsIOjTt], (/ 36 [ 0 , 27 r] 


d^{SpM^ gpM^^g) h dp dO 


n \ ClpkV{pk) 


= -u/ 3 'p{l cos pMid^Ii dp)v{p)^ ■ ^ {gM 2 M 2 ){ 9 , 9 k) 

Ae[o,2H . Cip 


k, k—h=:kl 


u(3' 

Hi 


Vip)^- Y1 {9M2M2){9,9k)cipkV{pk). 


k, k—h=±l 


(64) 
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with {gM2M2){9,0k) defined in (58) and where we nse the relation: 


/ cos dcp = — / sin if dcp 

J (pG[0,27r] J ip^[0,27r] 

= — sin^ ip Ml dip = — 

J ^G[ 0 , 27 v] 

obtained by integration by part and using (45). Using (61) and (64) and dividing by 

= I 


(/3S[0,27r] 


Sin 


pMiIidp) we get (56). 


Since I 2 does not depend on p, the contribntion of the last term to (60) for k = 2 


vanishes. From eqnations (24)-(25) and (34), we obtain the expression; 

TpM^^s{x,p,0,h,t) = - \^{p) X NpM-s{x,e,h,t) ■ z, 

NpM^^g{x,6,h,t) = p' sgn{k - h)j'^];^_ _{x,e,k,t). 

k, k—h=:kl 


Therefore, nsing (63) and since V{ip) = cos (</9 — p)V{p) + sin((y 9 — p)V{p)^ and Mi is 


even in we get after integration by parts: 


' 0G[O,7r], £/?G[0,27r] 


de{TpM- gpM^fi) h dp dd 


= p{ cospMidp)V{py 

J ip£[0,27r] 


with {gM2M2d0l2){9,9k) defined in 


p' Y {9M2M2del2) {9, h) ciPkV{pk) 

k,k—h=^l 


. ( 66 ) 


dividing by —C 4 = 2^2 ( 


6S[0,7r] 


sm 


29M2l2d9), we get (57). 


Using (62) and (65) into (60) for k = 


2 and 

□ 


4.5 Macroscopic equilibria 

One simple macroscopic eqnilibrinm consists in layers with the same vector velocity fields 
(or opposite vector field). In that case, inclinations have no impact on the dynamics and 
are simply transported by the velocity flow. In particniar, we have: 

Proposition 6. For any p e [0,27r], {ik)k& ^ {0,1}^ and {9k)k£Z ^ [0,7r]^, the homoge¬ 
neous functions 


Vh G Z, V(a:, t) G x M, p{x, h, t) = {—iy*^p, 9{x, h, t) = 9^, 

define a homogeneous macroscopic equilibria. 


The case of opposite vector flows may be non-stable since a small deviation from the 
eqnilibria leads to the global alignment of the layers. In particniar, nnmerical simnlations 
(see section 5.1.2) captnre only eqnilibria with the same velocity in each layers. The 


qnestion whether other stable macroscopic eqnilibria exists remains open. 
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5 Numerical experiments 


In this section, we compare nnmerical simulations of both the microscopic and macro¬ 
scopic models. The numerical methods are variations of those presented in Ref. 20 : the 


microscopic model is solved with an implicit scheme and the macroscopic model is solved 
using the splitting method. The two methods are detailed in 0 


5.1 Homogenous simulations 

5.1.1 Convergence to equilibria. 


We consider 3 layers and each layer contains 600 particles. The particle positions are 
uniformly randomly distributed on the square [0, L^] x [0, Ly] with = Ly = 1, the par¬ 
ticle velocity direction angles ip are uniformly randomly distributed on M/[0, 27r] and the 
particle inclinations 6 are uniformly randomly distributed on M/[0,7r]. We first consider 
an homogeneous test-case: the interaction radii Ri, R 2 and R^ all equal L^l2 and each 
particle thus interacts with (almost) all the other ones. 

We first consider non-interacting layers supposing h > 2R. In Figure we plot the 
distributions of ip and 6 of the three layers. We also plot the von Mises distributions: 






ie{l,2,3}. 


where ip^ and 0^ are the mean velocity angle and mean inclination angle of the particle of 
the £-th layer: 


= 


hj=i 




V p 

2^j,hj=e^ 


2iej(t) 


_ Q2i9j{t)\ 

5 ' 


Both velocity and inclination distribution are in good agreement with the von Mises 
distributions. In Figure we present the time evolution of the mean angles ip^ and 9^. 
Since there are no layer interactions, these mean angles are almost constant in time up 
to stochastic fluctuations. 


5.1.2 Interactions between layers. 

We still consider 3 layers but the number of particles per layer equals 25000. The particle 
positions are uniformly randomly distributed on the square [0, Lx] x [0, Ly] with Lx = Ly = 
1 and the interaction radii Ri, R 2 , R 3 equal 0.02: there are in average 31 neighboring 
particles. We include layer interactions: we suppose that the inter-layer distance h equals 
the particle radius R = 0.02. The layer interaction coefficients are chosen as follows: 
/i = 3 and u = 0.5. Particle velocity and inclination angles are randomly distributed 
according to their respective von Mises distribution. Initial mean velocity angle and 
mean inclination angle for the three layers are chosen as follows: 

= - 1 , = 1 , = - 2 , ( 66 ) 
9^ = 0, 9^ = -0.9, 9^ = -0.8. (67) 

We also perform a time rescaling in the microscopic model. Let £ > 0 and consider the 
following microscopic parameters: 

= z//e, D‘^ = D/e, = K/e, 5^ = 5/e, 


and (3'^ = e/i. 


( 68 ) 
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Velocity angle density, t =10.00 Orientation angle density, t =10.00 




Figure 2: (Homogeneous case, non-interacting layers) Left: Distributions of velocity angle 
ip at time t = 10. Right: Distributions of inclination angle 6 at time t = 10. Numerical 
parameters: u = 4, D = 0.3, K = 4, 6 = 0.3, At = 10“^. Number of particles: 600 per 
layer. 




Figure 3: (Homogeneous case, non-interacting layers) Left: Mean velocity angles as 
function of time. Right: Mean inclination angles r as function of time. Numerical 
parameters: u = 4, D = 0.3, K = 4, 6 = 0.3, At = 10“^. Number of particles: 600 per 
layer. 


with e: = 0.1. Consequently, we choose a macroscopic time scale. We compare particle 
simulations with macroscopic simulation. For the macroscopic model, we thus consider a 
constant initial density in each layer given by: 


(Number of particles) ^ forte {1,2,3} 

LxLy 


and the uniform initial values of and 9^ given by (66) - (67). The macroscopic layer 
interaction coefficients are given by: 


(3' = /?, jj,' = pRiiRl, 


with no e, since the particle simulation already consider the macroscopic time scale. 
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Fig, I (top) depicts the time evolution of the mean velocity and mean inclination 
for both the microscopic simulation (dashed line) and macroscopic simulation (continu¬ 
ous line). Microscopic simulations are averaged of 20 particle simulations. Let us hrst 
describe the macroscopic dynamics. The dynamics can be split into two steps: during 
the hrst step, up to time f ~ 2, layer 2 mainly interact with layer 1 since the overlap 
function is more important between this two layers. This leads to a hrst relaxation dy¬ 
namics that make the mean velocities of these two layers align. Then, during the second 
step, after time t = 2, interactions between layers 2 and 3 becomes predominant and a 
second relaxation dynamics occur that leads to alignment of the three layers. We then 
note that microscopic and macroscopic simulations coincide during the hrst 1.5 time unit 
(including the hrst relaxation mechanism). This conhrms that the macroscopic model 
captures the right interaction time scale. However, we see that, due to the hnite num¬ 
ber of particles, stochastic huctuations make the second relaxation occur earlier around 
time t = 3 (micro) instead oi t = 3.5 (macro). Looking at Figures]^ (bottom), where 
10 particle simulations are plotted and compared to macroscopic simulation, we see that 
the time of the second relaxation depends on the simulation and occurs always before the 


macroscopic relaxation. As noted in section 4A, once the particle velocities of the three 
layers are aligned, homogeneous inclination angles per layer dehne equilibria. Therefore, 
the time of the second relaxation step strongly determines the hnal inclinations. This 
explains the large deviation between macro and micro simulations after time t > 4 as 
regards inclinations. 

We now conserve the same parameters except that h = 0.0205 > R = 0.02 and, 
consequently, some inclination conhguration prevent layers from interacting. In Fig. 
we plot the time evolution of the mean velocity angle and mean inclination angle. We 
observe that, in the macroscopic simulation, a slight increase of the inter-layer distance 
results in large time translation of the second relaxation step going from t = 3.5 (Fig. 
1^ to f = 7.5 (Fig. [^. This highlights the meta-stability of the system between the two 
relaxation steps. Concerning the particle simulations (in dashed lines), the slight increase 
of the inter-layer distance does not result in a so much increase of the second relaxation 
time (it goes from t = 3 (Fig. to only t = 3.5 (Fig. [^). The second relaxation 
thus occurs two times earlier than predicted by the macroscopic simulation. Indeed, due 
to stochastic fluctuations, some particles interact instead of remaining in non-interacting 
conhguration. Therefore, this is the stochastic huctuations that impacts the long term 
dynamics of the model. 


5.2 Inhomogenous simulations 


We now consider 3 layers on the square domain [0, L^] x [0, Ly] with = Ly = 10. As 
in (^, we are interested in Taylor-Green vortex initial condition. Initial densities are taken 


uniform equal to 10'^. Velocity angles are given by 




u{xe,ye) 

vixe,ye) 


for i E {1,2,3}, 
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Figure 4: (Homogeneous case, interacting layers, h = R) Comparison of macroscopic 
(continuous line) and microscopic (dashed line) simulations. Left: Mean velocity angles 
as function of time. Right: Mean inclination angles r as function of time. Up: 
comparison with the average of 20 microscopic simulations. Down: comparison with 10 
microscopic simulations. Alignment interaction parameters: z/ = 4, D = 0.1,iC = 4, 
5 = 0.1, Ri = R 2 = R. Layer-interaction parameters: h = 0.02, R = 0.02, [3 = 0.5, 
R 3 = R. Microscopic parameters: 25000 particles per layer, /x = 3, At = 1 x 10“^. 
Macroscopic parameters: Ax = Ay = 0.5, At = 1 x 10“^. 


where w denotes the angle between vectors tc G and ( 1 , 0 )^, (u, v) is the vector dehned 
by: 

uix, t/) = ^ sin cos ^ sin cos ^ sin (^^x) cos (^^y^ , 

!,) = - i cos (Ji) sin (j!,) - t cos (APj i cos (^i) sin (^j/) . 

and {x£, yi) are translation of (x, y) : 

xe = x-tl mod y^ = y -ty mod Ly, 
with translation vectors given by: 

{tl tj) = (0, 0), {tl, tl) = (2,2), tl) = (5, 2). 
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Figure 5: (Homogeneous case, interacting layers, h > R) Comparison of macroscopic 
(continuous line) and microscopic (dashed line) simulations. Left: Mean velocity angles 
as function of time. Right: Mean inclination angles r as function of time. Up: 
comparison with the average of 20 microscopic simulations. Down: comparison with 10 
microscopic simulations. Alignment interaction parameters: z/ = 4, D = 0.1,iC = 4, 
(5 = 0.1, = i ?2 = R- Layer-interaction parameters: h = 0.0205, R = 0 . 02 , /3 = 0.5, 

R 3 = R. Microscopic parameters: 25000 particles per layer, /x = 3, At = 1 x 10“^. 
Macroscopic parameters: Ax = Ay = 0.5, At = 1 x 10“^. 


As regards the velocity initial condition, each layer is thus the translation of a normalized 
Taylor-Green vortex. Consequently, layers velocity helds are not initially the same and 
alignment dynamics should arise. Finally, inclination angles are taken uniform with the 
same value as the previous test-case (67). 


5.2.1 Non-interacting layers 

We hrst consider non-interacting layers : h > 2R. This test-case thus reduces to a 
simulation of the SOH model. In Figure we plot the space distribution of density, 
velocity and inclination for both the particle (left hgures) and macroscopic simulations 
(right hgures) for layer 2 (hrst and third layer are identical up to translation). For the 
particle simulation, we consider that each layer contains 10^ particles. The interaction 
radii i?i, R 2 and R 3 equal 0.04. Therefore, the particles have in averaged vrR^lO^ ~ 5 
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neighboring particles in each layer at the beginning of the simnlation. We consider the 
same rescaling (68) with e = 0.1. Densities are computed on a grid with space steps 


equal to Ax = Ay = 0.2 and as regards the particle simulation, they are averaged over 20 
runs of the test-case. We use the same time step for both macroscopic and microscopic 
simulations. 

In Fig. I^(top), we observe clustering for both microscopic and macroscopic simulation 
in region of negative divergence flow. The density in layer 2 has maximal value equal to 
||p 2 i|oo = 6520 for the microscopic simulations and ||p 2 ||oo = 7547 for the macroscopic 
simulation. Consequently the number of neighboring particles is multiplied by a factor 10 
in these regions. We also observe that the vortex are better conserved with the particle 
simulations. Concerning the velocity held, the microscopic velocity is obtained by dividing 
the local momentum by the local density. The velocity vectors should be of size ci but 
this is not the case in low density regions due to the small number of particles. This partly 
explains the observed differences between the microscopic and macroscopic simulations. 
However, there is a quite good overall agreement between the two simulations. Fig. 
(bottom) is represented the cosine (in absolute value) of the inclination angle: as layer 
interactions do not occur, it remains uniform. 


5.2.2 Interacting layers 

We then consider interacting layers when setting h = R. The layer interaction parameters 
are taken equal to: /3 = 2, /i = 20. The other parameters are the same as previously. The 
parameter /i is chosen in such a way that the interaction coefficient fiRnRlpo ~ 0.5 is 
of the same order as in section 5.1.2 We compare the average of 20 particle simulations 


with one macroscopic simulation on Fig. [^and[^ 

Fig. [^depicts the density and the velocity vector field for the three superposed layers. 
As in the previous test-case, we observe that the microscopic and macroscopic simulations 
are in good agreement. The velocity vector fields of the three layers are mostly aligned and 
consequently the density concentrations are localized almost in the same regions. This 
is particularly true for the macroscopic simulations (right). Concerning the microscopic 
simulations (left), we still observe some differences between the layers. 

In Fig. 1^ we represent the cosine of the inclination angle (in absolute value). The 
inclination for particle simulations is obtained by computing the local mean inclination 


angle with formula (10). Due to layer interactions, the inclination is no more uniform. 


Contrary to the velocity vector field, we observe large differences between the macroscopic 
and particle inclinations. This could be a consequence of the differences pointed out in 
section [5.1.2[ Note that regions with aligned inclinations do not necessarily match regions 


of uniform densities. Finally, the inter-layer interactions on inclinations affect in return 
the density and velocity vector field. This is particularly clear when comparing the density 
of Layer 2 with the one of the non-interacting test-case (Fig. (top)) for the macroscopic 
simulation. All the symmetries inherited from the initial distribution have been diluted 
by the inclination/velocity inter-layer interactions. 
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Figure 6: (Inhomogeneous case, non-interacting layers, velocity) Top: Velocity vector 
field. Down: Cosine of the inclination angle (in absolute value). Left: Particle simulation 
with £ = 0.1 (averaged on 20 simulations). Number of particles: 10^ per layer. At = 
1 X 10“^. Right: Macroscopic simulation. At = 1 x 10“^, Ax = Ay = 0.2. Alignment 
interaction parameters: u = 4, D = 0.1, K = 4, 6 = 0.1. Layer-interaction parameters: 
h = l> R = 0.02, R^ = R^ = R^ = R. 


6 Conclusion and discussion 

In this article, we have proposed an individual-based model of self-propelled disk-like 
particles interacting through alignment and volume exclusion. This model is intended 
to provide a framework for modeling collective sperm-cell dynamics. Particle motion is 
supposed conhned in two-dimensional planar layers. Particle interactions between nearby 
layers contribute to modify the disk inclinations, which generates a coupling between 
inclinations and motion. We have then derived a continuum model from this individual- 
based model. It describes the evolution of the local density, mean velocity direction and 
mean inclination of the disks in the various layers. Numerical simulations have shown a 
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good agreement between the continuum model and the individual-based one, but has also 
highlighted some differences. 

There are many possible directions to expand the current work and make it more 
realistic. At the individual-based level, the description of the agents and the interaction 
rules could be improved for a better account of actual sperm-cell motion. For instance, the 
assumption that all sperm-cells have the same constant velocity is obviously unrealistic. 
In any semen sample, there is always a certain proportion of dead sperm-cells and of 
less motile ones. This could be accounted for by allowing the particle velocities to span 
a certain range of values. The shape of the head could be improved from the current 
inhnitely thin disk to hnite thickness ellipsoids. The inclination interaction could involve 
a density dependency as it is more difficult to £t actual disks in one layer if they are 
inclined towards the plane than if they stand vertically. Finally, one could also imagine a 
process by which particles would change layers. At the level of the continuum model, one 
major improvement should be to add a random fluctuation term in order to account for 
hnite system size effects, similar to Ref. [^. We believe that adding such a term would 
help achieve a better match between the continuum model and the individual-based one. 
Other improvements would consist in adding a spatial diffusion to retain some of the 
nonlocality of the alignment interaction, similar to Ref. 11 or adding a layer-changing 


term. Finally, for both the individual-based and continuum models, the model parameters 
should be calibrated by close comparisons with biological data. 

A Nematic alignment 


Nematic alignment consists in alignment with the mean direction. The mean direction 
can be dehned as the eigenspace of the averaged projection matrix: 


cos 9 
sin 6* 


mdo, 


sm 


corresponding to the largest eigenvalue. The eigenvector (cos 0, sin 0)^ satisfies: 

^cos 

\ f I H \ riH I _ 1 V I 

ij 0 ^ ^ ^ 

Easy computations lead to the following relation: 


cos 6 

sin 9 j ^ \ sin 9 


= 0 


(69) 


(70) 


sm{2{9 - 9)) f{9) d9 = 0, 


(71) 


that can be written also as follows: 

'"cos 29 


I cos20 , „ 

sin20 / ^ lsin20 * “ 


The mean direction then satishes: 


cos 29 
sin 29 


cos 29 
sin 29 


cos 29 
sin 29 


(72) 


fdO 


fdO 


(73) 
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Then 6 is defined modulo n/2. Therefore, the orthogonal vectors (cos 0, sin 0)^ and 
(cos {6 + 7 r/ 2 ), sin(0 -|- 7 r/ 2 ))^ are both eigenvectors. Using the relation 



/cos 29\ 
\sm 29 J 


fd9 


„ /cos 29\ 
Vsin20 J ’ 


the averaged projection matrix writes: 


where C 



f cos 29\ 
\sin 29 J 


fd9 , 


1 A ± C cos 29 ±C sin 29 \ 

2 \ ± C sin 20 1 T C* cos 29J ’ 


(74) 


(75) 


whose eigenvalues are given by (1 ± C)/2. Consequently, the eigenvector (cos 0, sin 
corresponding to the largest eigenvalue involves the angle 9 satisfying: 


I' cos 29\ 
Vsin20 J 


r /cos 20 \ 

Jo Vsin2gy) 

r /cos 2 A 

Jo 


fd9 


fd9 


9 is now dehned modulo vr. We recover the dehnition of equation 10 


(76) 


B Coefficients of the macroscopic model 


From Ref. [^, coefficients ci and C 2 are positive and bounded. We have the following ex¬ 
pansion (in the limit ki —)■ 0 , that corresponds to large diffusion compared to alignment): 

Cl = -^ -L 0{kI), C2 = — -L 0{kI), 

The following proposition asserts similar results for coefficients C3 and C4. 

Proposition 7. Constants C 3 and C 4 can be written: 


C3 = 1 - 


vr 


-, C4 = 1 - 


vr 


^Kicosudu^ ^ —K,I cos udu^’’ ^K2COs2udu'^ ^ — K,2COs2udu'^ 


In particular, they are positive and lower than 1. We have the following Taylor expansions, as 

Ki — y 0 and Kj 2 — ^ 0 - 

C 3 = f + 0 (ft?), c, = ^ + 0{4). 

Proof. By integration by part: 


C3 = Ki / d^MJi dip = -Ki / Mid^Ii dip 

^V^G[0,27r] J ip^\fd.,2'R\ 


Then, using expression (51), we easily get the expected expression for C 3 . The same 


manipulations lead to the expression for C 4 . The positivity of C 3 and C 4 then results from 
the Cauchy-Schwarz inequality. □ 


Moreover the weight functions {gM 2 M 2 ) and {gM 2 M 2 dgl 2 ) are also bounded. 
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Proposition 8. We have: 

0 ^ {gM2M2){eA) ^ 1 , \{gM2M2del2){e,h)\ ^ 2 . 

We have the following Taylor expansions: 

r _ _ Jfj Jfjf 

{gM2M2)i9,ek)= / gie + e,e, + e') -+ 0 (^: 2 ), 

Je,e'e[o,iT] ^ ^ 

f - - d 6 dO^ 

{gM 2 M 2 d 0 l 2 ){O, 6 k) = —k ,2 / g {6 + 6 , 6 k + O') cos26 -I"0 (k2 )- 

J9,0'&[o,n] ^ ^ 

Proof. The inequalities for {gM 2 M 2 ) results from the bounds 0 ^ ^ 1. Using expression 


(52), we have : 

{gM2M2dgl2){6,6k) 

g {6 + 6 , 6 k + 6 ') ( — M 2 { 6 )+ 

71 


'9,e'e[o,n] 


= -{gM2M2){6,6k) 


^K,2COs2udu^ ^ — K.2COs2udu^ 


TT 


M2{6')d6d6' 


r _ _ 

(J- ^n,cos2udu^ (J- e-.,co.2udu^ Jgfi'elO,^] ^ ^ ^ \ ’ 

In the last expression, both terms have absolute value lower than 1. 


□ 


C Numerical schemes 

For the sake of completeness, we here recall the numerical scheme used for the numerical 
simulations. 


C.l Microscopic equations. 

The deterministic part of equations ([^ and ([^ can be written: 

dV{<pt) = Pvw,)^ [>'(r(^r) - '''(»>.))] dt, 

= (ie”"*) • (2A'(e““* - e”"*) + 2 Tt{ie“'‘^)) 

Following (Annex C), we then use the scheme: 

= g)k +2 V((p))) — B + \/2DAtek, modulo 27r, 


= 6 k + \26'l - arg(C')] -L VdAtSk, modulo vr. 


with C = + At 

'I 

where e^. and ef, are random variables with standard normal distribution, arg( 2 ;) G (0, 27r) 
denotes the argument of the complex number z G C and we recall that v G (0, 27r) denotes 
the angle between the vector n G and the vector (1,0)^. 
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C.2 Macroscopic equations. 


Multiplying equations (56) and (57) respectively by and we get: 


p{d,vm + cc^{y{p) ■ V.)V{p)) + Py^^y 


—v^p 

Ki 


= p 




u(3' 

C3 


{gM2M2){9,9k)c^pkV{yk) 


fc, k—h=^l 

p{dte^^^ + cc,{V{p)-V^)e^^^) = 
1 


C4 


(cipl^(<p))-^ ■ ^ sgn{k - h){gM2M2dgl2)i9,9k)cipkV{pk) 


2ie 


2ie 


(77) 


(78) 


/c, k—h=:kl 

To solve this system, we introduce a relaxation model 

dtp -L cci Vx ■ {pv) = 0, 

dtpv + CC 2 Vx ■ (pv <^v) -\ -VxP 

Ki 

= — 'y] {gM2M2){9{u),9{uk))ciPkVk--{l-\v\^)v, 

k,k-h=±l ' 

2 1 

dtpu + cciVa; ■ (pn ® u) = — (cipn)"*" ■ N -(1 — \u\Pu, 

C 4 p 

with 

N = p! ^ sgn(fc - h){gM2M2del2){9{u), 9{uk)) cipkVk, 

fc, k—h=^l 


( 79 ) 

(80) 
(81) 

(82) 


and where 9{u) G [0,27r[ denotes the angle of the vector u G In the limit p —)■ 0, 
the solution (p,v,u) to (79)-(80)-(81) formally converges to (p, 17(<p), (cos20,sin20)^), 
solution of equations (55) and (77)-(78). 

Equation (79)-(80)-(81) is numerically solved using a splitting method. We first solve 


the conservative part (with a Roe-like method [^), we then add the source term and we 
finally solve the relaxation part. For the last step, we just perform a renormalization of 
the vectors. This kind of scheme has been validated. In particular, it captures the correct 
discontinuous solutions of the macroscopic model corresponding to the solutions of the 


microscopic simulations. For more details, we refer to Ref. 20 


The characteristic velocities of the conservative system in the x direction are: 

7l = C (^C2Vx + \/ci/ki + vlc2{c2 - Ci)j , 72 = CC2Vx, 

= C {c2Vx - \/ci/Ki +n2c2(c2 - Ci)j , 74 = CCiU^,, 75 = CCiV^. 


73 


where denotes the first component of the vector v. As noticed in Ref. 10 , the conser¬ 
vative part of the equation is hyperbolic under the condition: 


Ci/ki +n^C2(c2 - Cl) > 0, 
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That is the case for the parameters chosen in section To ensnre the stability of the 
conservative step, the time and space steps At and Ax have to satisfy the CFL condition: 


max lyd At < Ax. 

l<i<5 

To ensnre the stability of the sonrce term step, we choose At small enongh to ensnre: 


At 

At 


ul 

C3 


{gM2M2){9{u),9{uk))cipkVk 


k, k—h=:kl 


< 2 , 


— (cipv)'^ ■ N 

C4 


< 2 , 


at each time step. 
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Figure 7: (Inhomogeneous case, interacting layers, velocity) Left: Velocity vector field for 
the particle simulation with e: = 0.1 (averaged on 20 simulations). Number of particles: 
10^ per layer. At = 1 x 10“^. Right: Velocity vector field for the macroscopic simulation. 
At = 1 X 10“^, Ax = Ay = 0.2. Alignment interaction parameters: u = 4, D = 0.1, 
K = 4, 5 = 0.1. Layer-interaction parameters: h = R = 0.02, /? = 2, /i = 20, Ri = R 2 = 
R 3 = R. 
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Figure 8: (Inhomogeneous case, interacting layers, inclination) Left: Cosine of the incli¬ 
nation angle (in absolute value) for the particle simulation with e: = 0.1 (averaged on 
20 simulations). Number of particles: 10® per layer. At = 1 x 10“^. Right: Cosine of 
the inclination angle (in absolute value) for the macroscopic simulation. At = 1 x 10“^, 
Ax = Ay = 0.2. Alignment interaction parameters: u = 4, D = 0.1, K = 4, 6 = 0.1. 
Layer-interaction parameters: h = R = 0.02, (3 = 2, y, = 20, Ri = R 2 = Rs = R. 

























